FERMILAB-PUB-98/181-T 
June, 1998 



An Efficient Algorithm for QCD 
with Light Dynamical Quarks 

A. Duncan^, E. Eichten^, and H. Thacker^ 

"'^Dept. of Physics and Astronomy, Univ. of Pittsburgh, Pittsburgh, PA 15260 
^Fermilab, P.O. Box 500, Batavia, IL 60510 
^Dept.of Physics, University of Virginia, CharlottesviUe, VA 22901 



Abstract 

A new approach to the inclusion of virtual quark effects in lattice QCD 
simulations is presented. Infrared modes which build in the chiral physics in 
the light quark mass limit are included exactly and in a gauge invariant way. At 
fixed physical volume the number of relevant infrared modes does not increase 
as the continuum limit is approached. The acceptance of our procedure does 
not decrease substantially in the limit of small quark masses. Two alternative 
approaches are discussed for including systematically the remaining ultraviolet 
modes. In particular, we present evidence that these modes are accurately 
described by an effective action involving only small Wilson loops. 
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1 Introduction 

Recent studies of the origin and role of exceptional configurations [|], ^ leading to ex- 
tremely noisy hadron correlators in quenched lattice QCD for light quark masses have under- 
lined the importance of nonlocal topological fluctuations in determining the chiral physics of 
the quenched theory. It has been known for a long time Q| that such fluctuations completely 
alter the behavior of the theory in the light quark mass limit once the quark determinant is 
included in a full dynamical calculation. Moreover, the appearance of spurious real modes 
of the Wilson-Dirac operator on finite lattices in association with such fluctuations ac- 
counts for the increasing frequency of exceptional configurations at stronger coupling and 
low quark mass. These connections all suggest that reliable calculations in the chiral limit 
of lattice QCD require an accurate treatment of the low cigcnmodcs of the quark Dirac 
operator. 

Singularities of quark propagators in the quenched theory are automatically regulated by 
corresponding zeroes of the quark determinant. However, this regularization is only effective 
if the low eigenmodes of the quark Dirac operator are treated precisely: in particular, 
valence and sea quark masses must be identical. On the other hand, the high eigenmodes 
(corresponding to imaginary masses for scales well above the QCD scale up to the lattice 
cutoff) when integrated out contribute to an effective gauge-invariant gluonic action which, 
for physics on small momentum scales, simply amounts to a redefinition of the scale of the 
theory. To the extent that the ground state hadron spectrum involves hadronic bound states 
with constituent quarks all off-shell on the order of Aqcd, it therefore seems likely that high 
eigenmodes are simply irrelevant for spectrum calculations, even in full QCD. By writing 
the fermion determinant in terms of the hermitian operator H = 75 (7^ (A) — m) |^ we are 
able to deal with a completely real spectrum. Moreover, the individual eigenvalues have a 
direct physical interpretation as a gauge-invariant measure of off-shellness of the quark fields 
(to see this, we recall that for the free continuum theory, the eigenvalues of H are simply 
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±-\/p^ + for a quark mode of Euclidean momentum p). 

We shall argue in this paper that a separation of low and high eigenmodes can in fact be 
carried out in a practical way in unquenched lattice QCD calculations, leading to an efficient 
way of building in the important physics of the quark determinant in the chiral limit. Such 
a separation also corresponds to a completely gauge-invariant and smooth interpolation 
between the quenched and full dynamics of the theory. The procedure we propose also 
yields as a byproduct very detailed and useful information about the infrared spectrum of 
the Dirac operator which is known to be intimately related to the chiral physics of the theory 

and central to the overlap formulation ||^ of lattice QCD. 

In Section 2 we describe a regularized version of the fermionic determinant which in- 
terpolates smoothly between the quenched and full theory, in a way which allows for the 
selective inclusion of fermionic modes in a predetermined momentum range (typically from 
zero up to a given cutoff /i). This regularization is amenable to an analytic perturbative 
calculation in which the role of the high eigenmodes contributing to the full fermionic deter- 
minant can be clearly isolated. Such a regularization can be studied analytically in abelian 
4D gauge theory, where the /^-dependence of the determinant for large /i and low momentum 
is seen to reduce to a shift of bare coupling (or in the lattice context, of scale). 

In Section 3 we describe the results of some simulations in 2 dimensional lattice QED 
(QED2), which has proven to be an extremely useful testbed for exploring features of the 
Dirac- Wilson spectrum in lattice gauge theory. Here, and henceforth in all the numerical 
simulations, we employ a regularization of the fermionic determinant in terms of a sharp 
mode cutoff which is physically equivalent to the smooth regularization of Section 2 but 
suitable for numerical implementation in large systems. It is shown that the fluctuations 
of the full fermionic determinant in an exact dynamical simulation of QED2 are essentially 
restricted to a small fraction (for the lattices studied here, essentially the lowest few per- 
cent) of the spectrum. Comparisons of pseudoscalar correlators computed in the quenched 
and full dynamical theory are made with an approximate simulation where only the low- 
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est ten percent of the eigenvalues of the hcrmitian operator ^^{Ip— m) are included in the 
fermionic determinant. The truncated determinant simulations essentially reproduce the 
full dynamical results. Two characteristic features of unquenched gauge theory, the sup- 
pression of topologically nontrivial sectors and the breaking of the string due to shielding, 
are also illustrated using the truncated determinant approach in QED2. Of course, in the 
case of QED2, the superrenormalizability of the theory implies that the high cigenmodes 
are basically inert, in distinction to the case in 4D gauge theory where these modes will 
necessarily introduce a further logarithmic rescaling due to the variable screening effect of 
virtual quark/antiquark pairs at different length scales. 

In Section 4 we describe in detail the algorithm we have employed for the simulations 
of full QCD (on a 12^x24 lattice at P—5.9 and inverse lattice spacing a^^ =1.78 Gev for 
the quenched theory) with a truncated determinant. The Lanczos procedure allows reliable 
extraction of Dirac cigenmodes up to energies ^370 MeV, certainly enough to include the 
essential low energy chiral physics of QCD. Moreover, the Lanczos procedure extracts the 
needed small eigenvalues rapidly as the spectrum is relatively sparse there. Unlike the case 
of propagator inversion, the Lanczos method is stable even in the presence of very small 
eigenvalues provided these are not too dense. We also discuss some aspects of the Monte 
Carlo dynamics (acceptance rate and equilibration time) for our update procedure, in which 
pure gauge heat bath sweeps alternate with Metropolis accept/reject steps for the truncated 
determinant. A crucial point is that we do not see a dramatic fall in the acceptance rate of 
our procedure as we go to lighter quark masses. 

In Section 5 we present the results of our truncated determinant simulations of QCD4. 
The initial study involves runs on a 12*^x24 lattice at f3—5.9 and at three kappa values 
(0.1570, 0.1587 and 0.1597), reaching in the lightest case a pion mass on the order of 
280 MeV. Pseudoscalar meson masses are measured and a value for the critical hopping 
parameter extracted. The inclusion of 100 quark cigenmodes (all modes up to ~370 MeV) 
eliminates the necessity for considering quenched chiral logs in the chiral extrapolations. 



4 



The topological charge distribution is measured for different quark masses and compared 
with the quenched result. As expected, nonzero topological charge is strongly suppressed in 
the light quark limit. Measurements of the string tension reveal clearly a screening of the 
quark-antiquark potential from the virtual sea quarks, although the lattice used is still too 
small to allow us to see the asymptotic flattening expected at large distances. 

In Section 6 we show that the high momentum modes can be included in a precise way 
by a combination of the truncated determinant and multiboson methods |^ . The procedure 
suggested in this paper is in a sense exactly complementary to the multiboson approach of 
Liischer. The latter approach treats the high eigenmodes of the Dirac operator very well, 
but necessarily introduces errors whenever small eigenvalues are present. In the chiral limit 
such modes become frequent and in fact dominate the chiral physics. Here we propose 
treating these modes as precisely as possible. Another approach to the inclusion of the high 
modes, a loop Ansatz for the short distance piece of the quark determinant, is also discussed 
in this final section. Such an Ansatz, involving relatively short Wilson loops (up to length 
6), is shown to give a very accurate description of the high end of the quark determinant. 
Finally, in Section 7 we summarize our conclusions. 

2 Truncated determinants in gauge theory 

The separation of low and high eigenmodes in the fermionic determinant can be accomplished 
in an analytically convenient way by smoothly switching off the higher eigenvalues above a 
sliding momentum scale /i. Given a matrix M then det(tanh(Al//i)) reduces to unity for /i 
much below the smallest eigenvalue of A4 while reproducing the full determinant (up to an 
irrelevant multiplicative factor) for fi much above the highest eigenvalue. For a gauge theory, 
defining the hermitian operator H = j5{I^{A) — to), then the effective action obtained from 
integrating out each flavor of fermion of mass m can be regularized as the logarithm of the 
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smoothly truncated determinant 

V{^,) = itrln(tanh(:^)) (1) 

= tr{ln(l-e-2«'/''')-ln(l + e-2ff'/^')} (2) 

This definition allows an analytic calculation of the regularized determinant in weak 
coupling perturbation theory, which describes the /^-dependence of ©(/i) for /x well above 
the QCD scale. The calculation can be carried out for a nonabelian lattice regularized 
theory, but wc shall illustrate the procedure here for the case of a continuum 4-dimensional 
abelian gauge theory. Note that 



H'^ = K0 + K1+K2 (3) 

Ko = D+m^ (4) 

K, ^ {-i^,^} (5) 

K2 = A^A^ (6) 



To second order in weak coupling perturbation theory, we may compute ©(/u) by expanding 
to first order in K2 and to second order in Ki. The calculation is lengthy but straightforward 
(details will be given elsewhere)- here we quote the result only. Expressed in terms of 
momentum space fields, one finds 

© = 1 (04/?(^'' ^' ^^)A^ik){k^5^, - k^K)A,{-k) (7) 

The contribution of the high modes can be studied by taking /x large compared with the 
quark mass m and, in the nonabelian case, with the QCD scale. Then these modes affect 
the low energy physics (i.e for fc << /i) through the low momentum limit of Explicit 
calculation gives 

(3{e,m,n)c,-l-H^) + 0{k') (8) 

which exactly corresponds to the expected /x-dependence of the screening shift in the running 
coupling induced by virtual fermionic modes in the momentum range up to /i. 
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The decoupling of the high fermionic modes suggests that lattice QCD calculations 
performed at weak enough coupling should be insensitive to the fluctuations induced by 
eigenvalues of the Dirac operator much above the QCD scale, except for an overall shift in 
the scale of the theory induced by renormalizations of the coefficients of the low dimension 
operators making up the effective pure gauge action. In particular, dimensionless ratios of 
physical quantities should fairly soon become insensitive to inclusion of higher modes in the 
fermionic determinant. In a superrenormalizable theory like QED2, this inscnsitivity should 
even be apparent in dimcnsionful quantities, as we do not have a logarithmic running of 
scale in this case. 

3 Truncated Determinant Algorithm in QED2 

Abelian gauge theory in 2 space-time dimensions (the massive Schwinger model) has proven 
to be a marvellously manageable testbed for exploring in detail |^ the spectral properties 
of the Dirac- Wilson operator. The computational expense of performing even exact update 
full dynamical simulations is relatively slight, essentially full information on the spectrum 
can be obtained configuration by configuration, and the system mimics, at least qualitatively, 
many of the topological and chiral properties of 4 dimensional QCD. This model also turns 
out to be a very useful starting point for investigating the relative importance of the infrared 
and ultraviolet ends of the Dirac spectrum in a full dynamical lattice simulation. 

Although the calculation of all the eigenvalues of and hence of as defined in 

the previous section, is perfectly feasible for 2D QED, the restriction of practical numerical 
techniques for the much larger matrices of 4D QCD to the low-lying eigenvalues suggest the 
use of a simpler truncation of the determinant, in which the lowest (in absolute magnitude) 
N\ positive and negative eigenvalues of H are included and all higher modes dropped. As 
we shall see in Section 5, precisely such a truncation scheme matches on exactly to a very 
accurate representation of the high end of the determinant in terms of an effective loop 
action. Labelling positive eigenvalues of H as rjn and negative eigenvalues as Cn (where the 
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index runs in the direction of increasing absolute magnitude) we define 

1 

V{N,) ^ ^J2^n{rjlO (9) 

n=l 
1 ^ 

where 2D is the dimensionahty of the discrete Wilson-Dirac matrix for the lattice theory. 
An exact full dynamical simulation would include the full (log) determinant ^{Nx) +'D{N\) 
in the effective gauge action. The extent to which the low eigenvalues determine the physics 
of the unquenched theory can be examined by comparing the fluctuations- in a dynam- 
ical simulation- of 'D{N\) with those of T>{N\) for various choices of N\ << D. These 
fluctuations are shown graphically in Fig. 1 , for 40 configurations generated in a full dynam- 
ical simulation using an exact update algorithm. In Fig. 2 the fluctuations are shown for 
40 configurations in a quenched simulation. The lattice used was 10x10 at /3=4.5 with a 
bare quark mass of 0.095. Evidently the fluctuations are essentially all confined to the low 
end of the spectrum. In the quenched case the size of the fluctuations at the infrared end 
is considerably larger than for the dynamical configurations, as configurations with small 
eigenvalues are suppressed once the determinant factor is included in the update procedure. 
The appearance of such configurations is intimately related to the exceptional configurations 
encountered in quenched calculations at strong coupling and/or small quark mass. 

This behavior suggests an approximate unquenched algorithm in which only 'D{N\) is 
used in the determinant part of the effective lattice action. If iV^ is chosen large enough, all 
the nonperturbative infrared physics will be properly included. Exceptional configurations, 
in which there is an anomalously low eigenmode of the Dirac- Wilson operator, are tamed 
in the expected way [ pT| , and the convergence of the procedure can be examined simply 
by repeating the run for increasingly large N\. The update algorithm we have chosen is 
very simple: a number (typically 5) of conventional Monte Carlo (Metropolis) sweeps are 
performed to obtain a new gauge configuration, a new value for 'D{N\) is calculated (for 
QED2 on a 10x10 lattice, we can easily obtain all the eigenvalues by direct diagonalization) 
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Figure 1: Fluctuations in T>{Nx),T>(Nx), dynamical configurations. For ease of 
visibility, the curves have been shifted vertically. 
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Figure 2: Fluctuations in 'D{Nx),'D{Nx), quenched configurations 



and compared with that for the preceding configuration. Then the truncated determinant 
factor is used to provide a Metropohs accept /reject criterion for the gauge configuration 
update. With N\—10, the acceptance ratio was typicaUy in the range of 50-75%. Measured 
quantities such as the pseudoscalar correlator decorrelated after a few configuration updates 
(the statistical errors shown include autocorrelation times computed from the data). The 
results of such a procedure for the pseudoscalar correlator are shown in Fig 3 for a quark 
mass (lattice units) of 0.05. For such a small mass, using Wilson fermions, the quenched 
functional integral is dominated by real pole contributions which appear in the simulation 
as wildly noisy values for the correlators (one finds in a typical run of 800 sweeps values 
for the pion propagator at zero time separation ranging from several thousand to zero, for 
a quantity averaging to order unity in the dynamical theory). Instead we have plotted the 
quenched results regularized by the pole shifting (or "modified quenched approximation" ) 
procedure of ||]. Evidently the full dynamical result for the pseudoscalar correlator is 
reached with only a small fraction (in this case, about 10%) of the eigenmodes included 
in the fermionic determinant. When only 5% of the eigenvalues are included, the result is 
intermediate between the quenched and full dynamical values. 

As we are working at a very small value of quark mass, it is also of interest to study 
the effect of the truncated determinant factor on the topological charge distribution and the 
string tension of the theory. For the quenched simulations on a 10x10 lattice at /3 =4.5, the 
topological charge, defined as 



(where 0p is the plaquette angle for plaquette P), is found to be concentrated at roughly 
integer values, with charges and 1 dominating. The histogram of topological charge values 
obtained from 800 quenched configurations is shown in Fig. 4. As low eigenvalues are 
introduced via the truncated determinant, the nonzero topological charge configurations 
are suppressed. Again, with N\~5, the resulting distribution is hardly distinguishable from 
the full dynamical result. 
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Figure 3: Comparison of quenched, full dynamical and truncated determinant sim- 
ulations: 10x10 lattice, 13=4.5, m=0.05 

The quark-antiquark potential determined for two different sea quark masses (bare mass 
0.06 and 0.10) is shown in Fig 5. The calculation was done in the truncated theory on a 
16x16 lattice at /3~4.5 using A^a=10 eigenvalues. Also shown is the exact result for the 
quenched theory and exact update full dynamical theory (at sea quark mass 0.06). The 
small eigenvalues properly represent the large loops induced by the determinant, leading to 
a breaking of the string at longer distances for light quarks. 

4 Calculating Truncated Determinants in Large Sparse 
Systems 



There are a number of techniques available for the extraction of a limited number of low 
eigenvalues of a large sparse linear system: most popular are conjugate gradient approaches 
[E2|, or the Lanczos technique p3|, suitably modified to guard for the appearance of spuri- 



ous eigenvalues |14 . The latter approach has been studied extensively by Kalkreuter |15|1 , 
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Figure 4: Topological charge distribution for quenched, full dynamical and truncated 
determinant simulations: 10x10 lattice, /3=4.5, m=0.05 




Figure 5: Quark-antiquark potential in QED2. 
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and has proven to be most well suited for the task at hand, namely, the accurate extrac- 
tion of a complete set of low- lying eigenvalues of 75 (i^ — m) up to an energy scale which 
ensures that all the important soft chiral dynamics of full QCD is included in the Monte 
Carlo simulation. Typically, on lattices of physically interesting volume in 4D QCD, this 
requires the determination of something on the order of 100 eigenvalues. One advantage 
of the Lanczos approach is that it can be pushed through to the determination of as many 
eigenvalues as desired, with a numerical effort which grows empirically as roughly the square 
of the desired energy cutoff (in the portion of the spectrum corresponding to the physical 
branch). In particular, on small lattices, it is relatively straightforward to determine the 
entire spectrum, which is useful both for diagnostic purposes and in studying the systematic 
effects of an algorithm based on a truncated determinant. 

The Lanczos technique is a standard part of the literature in Numerical Analysis (see, 
for example, jl^ ) so we shall give only a very brief review of the procedure here. Given a 
hermitian matrix H (in our case, this is just the matrix 75 (-/j/' — m) discussed previously), 
a series of orthonormal vectors vi,V2,V3... are generated from a starting vector wq = vi by 
the following recursion: 

Vn+l = Wn/I3„ (12) 

n 71+1 (13) 

a„ = {Vn.HVn) (14) 

Wn = {H - a,J)Vn - /3n-lVn-l (15) 

f3„ = \J {Wn,Wn) (16) 

where a„, /3„ are real numbers (by virtue of the hermiticity of H), and the initial conditions 
are (3o — l,vo — 0. It is straightforward to verify that the matrix of H, in the basis spanned 
by the Lanczos vectors Vn is tridiagonal, with the numbers ai, a2, .. down the main diagonal, 
and the numbers Pi, (32, ■■ on the first super- (and sub-) diagonals. If the Lanczos recursion 
is carried to order N, one is thus led to a NxN truncation of the original linear system, 
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and the eigenvalues of the resulting tridiagonal matrix T^^-* represent increasingly accurate 
approximants to the eigenvalues of the full system, provided only that the starting vector 
wq is not entirely contained in an invariant subspace of H. Typically, wq is chosen randomly 
and one expects that such a random vector will overlap nontrivially with all the eigenvectors 
of iJ. 

There are several features of the Lanczos procedure which appear at first sight problem- 
atic but which nevertheless turn out not to compromise its efficacy in the present application. 
First, degenerate eigenvalues of the original matrix H arc not properly handled (although 
methods have been devised for circumventing this drawback p3|). This turns out to be 
irrelevant in our QCD application as the generic spectrum of H = 75 m) for a typical 
gauge configuration encountered in the course of a Monte Carlo simulation is entirely non- 
degenerate. Secondly, the effects of roundoff error in the algorithm can be quite severe, and 
lead to the appearance of spurious eigenvalues and to the false duplication of real eigenval- 
ues. Fortunately a simple and effective cure for this problem, first suggested by CuUum and 
Willoughby proves to be practical in the gauge theory case [p_5l. However, it remains 
an unfortunate feature of the algorithm that the number of accurate eigenvalues extracted 
at level N of the recursion is typically considerably smaller than N. For example, on a 
12^x24 lattice at f3—5.9, the extraction of the lowest 100 eigenvalues of the Dirac operator 
typically requires on the order of 10,000 Lanczos sweeps. (The computational cost of a 
single Lanczos sweep is essentially that of the single multiplication incurred in produc- 
ing the next Lanczos vector Wn — Hv„ + ...) Finally, although the diagonalization of a 
tridiagonal matrix is conceptually trivial and efficiently implementable by scalar algorithms 
(e.g. by a QL implicit shift algorithm |l^), the parallel implementation of this procedure 
is not entirely trivial. In our simulations, this is essential to avoid a serious bottleneck in 
the simulation when tridiagonal matrices of order up to several tens of thousands must be 
efficiently processed. We describe below an elegant parallel approach to the extraction of 
the spectrum of T'^^K 
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In the case of the Dhac operator in QCD, the choice of the random vector wq used to 
start the recursion appears to be fairly innocuous (a local source seems perfectly adequate, 
for example) . An important check that the spurious eigenvalues are correctly identified 
and that the remaining "good" eigenvalues are sufficiently converged relies on the gauge- 
invariance of the individual eigenvalues of H which can be verified explicitly by recomputing 
the eigenvalues with varying degrees of gauge-fixing. We have performed extensive checks 
to ensure that the eigenvalue spectrum, and a- fortiori the truncated determinant 'D(Nx), 
is invariant (typically to at least 8 significant figures) under gauge transformations of the 
input configuration. 

The procedure we use to isolate converged eigenvalues of H involves two stages. First, 
spurious eigenvalues are identified by the procedure of CuUum and Willoughby- namely, 
eigenvalues of the tridiagonal matrix T^^'> are compared with those of the matrix r(^) ob- 
tained by deleting the first row and column of T^^'' and removing all common eigenvalues 
of the two systems. Secondly, converged eigenvalues are identified by requiring either du- 
plication of the remaining good eigenvalues or stability within a preassigned precision level 
when eigenvalues are compared at recursion level N — iVgap and N. Typically we insist on a 
precision level of at least I0~^ and choose A^gap=f 00. The above procedure requires the res- 
olution of the central part of the eigenvalue spectrum for four large tridiagonal matrices (of 
dimension N,N-1,N — A'gap, and N — A'gap-f , respectively). It is therefore highly desirable 
to perform these diagonalizations in a way that allows parallelization. 

The key to extracting the central part of the spectrum of a tridiagonal matrix in a parallel 
machine lies in the Sturm Sequence property of such matrices [^3|, valid provided none of 
the subdiagonal entries /3„ vanish, as is certainly the case for generic gauge configurations. 
Let Pn(A) be the secular determinant det{Tn^^ ~ XI) of the nxn principal submatrix xji^'^ 
of T'^^\ Then the number of eigenvalues of T^^^^ less than any preassigned A equals the 
number of sign changes in the sequence po(A),pi(A), ...pAr(A). As the Pn(A) can readily be 
calculated by a simple two term recursion relation, an obvious bisection procedure can be 
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Figure 6: Spectral density,/?=5.9,«;=0.1587 

used to determine tlie n'tli eigenvalue of T^^^ at any desired level of precision. Moreover, 
the task of extracting different eigenvalues can be assigned completely independently to 
separate processors, provided global access to the elements (ai,/3i) of T^^^ is arranged. 

In Fig 6 we show the spectrum of H for a typical configuration on a 12^^x24 lattice 
at /3=5.9, k;=0.1587. In this case an extended Lanczos recursion was carried out up to 
order A/'=78000, yielding 1478 converged eigenvalues in the central region of the spectrum. 
Converting the gauge-invariant eigenvalues of -ff to a physical energy scale (using the scale 
o~^=1.78 GeV from the charmonium spectrum), this corresponds to all quark eigenmodes 
up to 970 MeV. The energy reach as a function of number of Lanczos sweeps for this lattice 
is shown in Fig 7. In the simulations reported below, we have typically used 9500-12000 
Lanczos sweeps (with slightly different tunings of the CuUum-Willoughby procedure) and 
included the lowest 100 (i.e. 50 positive and 50 negative) eigenvalues of H in the update 
procedure. This cutoff corresponds to inclusion of quark eigenmodes up to an energy of 
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Figure 7: Lanczos convergence,/3=5.9,K;=0.1587 

approximately 370 MeV. Lanczos recursion to order N ~10000 requires (for a 12^x24 lattice) 
about an hour on 64 nodes of the Fermilab ACPMAPS system. 

The relaxation of the determinant from its typical quenched value to the equilibrium 
value appropriate for the system simulated with the truncated determinant is shown in 
Figs. 8 and 9. On a small lattice, (6^ at /3 = 5.7, fi;=0.1685) a dynamical run including 
N\ =30 eigenvalues (or up to about 500 MeV in quark mode energy) was performed, with 
a determinant update accept/reject every 3 heat-bath sweeps of the lattice. The resulting 
evolution of ViNx) starting from a quenched configuration is shown in Fig 8. The same 
evolution is shown for a run on a 12^x24 lattice at /3=5.9, k=0.1587 with N\ =100 in Fig 
9, and with determinant accept/reject every 2 heat-bath sweeps. 

The procedure we have used does not show a very strong dependence of the acceptance 
rate on the quark mass (fortunately!). For the heaviest quark mass studied at /? =5.9, 
K =0.1570 (slightly lighter than the strange quark), the acceptance rate for a run with 
determinant accept/reject performed every two gauge heat bath sweeps was 40%. For the 
lightest two quark masses (k=0.1587 and 0.1597) two separate runs were performed with 
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Figure 9: Determinant relaxation,/3=5.9,K=0.1587,A^A=100 
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determinant accept /reject steps separated by either one or two heat-bath sweeps. For the 
heavier case, k=0.1587 (corresponding to a pion mass of around 400 MeV), the acceptance 
was on the order of 37% for new configurations separated by 2 heat-bath sweeps, and 57% 
for configurations separated by a single heat-bath sweep. For the fighter mass, k=0.1597, 
(corresponding to a pion mass of around 280 MeV) the acceptance was 30% for new config- 
urations separated by 2 sweeps, and 43% for those separated by a single heat-bath sweep. 

5 Simulations in QCD4 

To test the efficacy of the truncated determinant approach to unquenched QCD we have 
performed some preliminary runs on a 12x^x24 lattice at P—5.9, for two degenerate flavors 
of dynamical quarks at hopping parameter values k = 0.1570, 0.1587 and 0.1597. For the 
heaviest mass, n =0.1570, propagators were computed for every fifth configuration (60 in 
all), and determinant accept/reject performed after every two heat-bath gauge updates. 
For the lighter two masses, we also performed parallel runs with determinant accept/reject 
after every heat-bath sweep and with propagators measured every tenth configuration. Our 
results are based on 104 propagators for k =0.1587 and 88 propagators for k =0.1597. As 
fluctuations in the large distance behavior at light quark masses grow, the run at k =0.1597 
is continuing and results with much higher statistics will be presented in a later work. 

Pseudoscalar meson masses were determined by measuring meson correlators with a 
smeared source and local sink and doing a fully correlated Euclidean time fit in the time 
window 8-11. At kappa values k=0.1570, 0.1587, and 1597, the pion masses were found to be 
(lattice units) = 0.339±0.011, 0.241±0.019 and 0.198± 0.069. The extrapolation to zero 
pion mass gives a critical kappa value of Kc w0.1602, slightly higher than the pure quenched 
value of Kc=0. 15972 found in previous calculations |jl^. With this critical kappa, the pion 
mass corresponding to our lightest case should be about 0.157 (lattice units), which converts 
to about 280 MeV using the scale determined for the quenched theory at this beta. We do 
not expect the scale in the truncated determinant simulation to be much changed from the 
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quenched theory as the shift in beta in fuU dynamical simulations (see [Q) is due to a large 
logarithm accumulated from quark modes all the way up to the lattice cutoff, whereas the 
truncated determinant included here only takes into account the infrared modes up to 370 
MeV. Higher statistics are presently being accumulated at the lightest quark mass, where 
the fluctuations in the large time meson correlators are largest. 

An important aspect of the low energy chiral dynamics of QCD is the response of the 
topological structure of the theory to the presence of light dynamical quarks. Integrating 
the vacuum expectation value of the U(l) axial anomaly in QCD yields immediately, in the 
continuum, the anomalous chiral Ward identity 



We may therefore define a topological charge on the lattice by simply evaluating the left- 
hand-side of the above equation configuration by configuration. The advantage of this 
definition is that the required information is already immediately accessible: the Euclidean 



H = "/^{I^ — m), the low eigenvalues of which were extracted in the course of the simulation. 
Namely, we have the following lattice definition of Qtop 



where the matrix dimension of H is N and Xi are the eigenvalues of H. 

As the larger eigenvalues occur roughly as equal and opposite pairs, this sum actually 
saturates quickly at the low end, and the 100 eigenvalues computed already in the course of 
the simulation suffices to determine Qtop to a few percent. An example of the convergence 
of this spectral sum, with the mode eigenvalues converted to a physical energy scale (recall 
that the inclusion of 100 eigenvalues corresponds to modes up to about 370 MeV), is shown 
in Fig |l^ for a typical configuration. 

With the definition (18), the qualitative effect of the quark determinant on the topo- 
logical charge distribution can readily be studied. Two effects are clearly visible in our 




(17) 



vacuum expectation value J 



X < ip{x)^^'4>{x) > reduces to the trace of the inverse of 




(18) 
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Figure 10: Convergence of the topological charge spectral sum for a typical config- 
uration 

data: 

(1) The real mode artifacts characteristic of quenched Wilson gauge theory, and which 
underly the increasingly frequent appearance of exceptional configurations as one goes to 
lighter quark masses, do not occur. Configurations with very small eigenvalues of H are 
suppressed by the determinant factor which is driven strongly negative for any such con- 
figuration. In Fig |l^ this effect is seen comparing the topological charge distribution for a 
dynamical run a,t n — 0.1587 (including the lowest 100 modes) with a completely quenched 
run at the same valence quark mass, Fig|l^. The distributions are broadly similar (although 
the dynamical one is slightly narrower) but the outlying points corresponding to very large 
values of Qtop (i-e to the appearance of an exactly real mode of the Wilson-Dirac operator 
very close to the chosen kappa value) are eliminated in the dynamical run. In fact in the 
quenched case there are several outlying points (the furthest out at Qtop = -91.6) not shown 
on the figure. For the dynamical run only charges IQtopI < 5 are seen. 

(2) Nonzero topological charge must be suppressed in the chiral limit of vanishing quark 
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Figure 11: The topological charge frequency distribution computed for k = .1587 
using 100 decorrelated quenched QCD configurations on a 12^x24 lattice at (3 = 5.9. 
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Figure 12: The same topological charge frequency distribution for 300 configurations 
generated by the truncated determinant algorithm with sea quark mass k = .1587 
on a 12^x24 lattice at (3 = 5.9. 
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Figure 13: The topological charge frequency distribution with sea quark mass k 
.1570. (Other parameters as in the previous figure.) 
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Figure 14: The same topological charge frequency distribution with sea quark mass 
K = .1597. The expected distribution from chiral perturbation theory is shown with 
diamonds. 
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mass, so we expect that the histogram of measured topological charges will narrow as one 
approaches Kc- This effect is shown in Figs. |l^ and |lj where the topological charge distri- 
bution is compared for two dynamical runs at the lowest and highest quark masses studied 
(k= 0.1597 and 0.1570). The narrowing of the distribution for the lighter mass is immedi- 
ately apparent. For comparison we plot the analytic result predicted by the chiral analysis 
of Leutwyler and Smilga in the light quark limit for two degenerate flavors: 

P(Q) = Iq {xf - Iq+, (x)/q_i (x) (19) 

with X = \VFlMl. Here V,F^,M^ denote the lattice space-time volume, the pseu- 
doscalar decay constant and the pion mass respectively, all in lattice units. For k=0.1597 we 
have taken M7r=0.15 and F^=0.07 (the latter number is extrapolated from high statistics 
quenched runs for this lattice [|l6j). 

The presence of low-momentum virtual sea-quark modes in the simulation should result 
in screening of the quark-antiquark potential extracted from Wilson loops at large distance. 
In Fig |l5|the potential obtained in the quenched theory on a 12^^x24 lattice at (3 =5.9 (200 
configurations) is compared with that calculated from our dynamical configurations at the 
lightest sea quark mass (k =0.1597). The effect of screening is clear although asymptotic 
flattening of the potential on this lattice occurs at distances where statistical fluctuations 
as well as finite volume effects dominate. 

6 Matching the High Eigenvalues 

Because QCD in four dimensions is only renormalizable, not super-renormalizable, the fluc- 
tuations of the ferniion determinant are significant at all physical scales. Therefore, unlike 
QED in two dimensions, it can not in general be sufficient to compute accurately only the 
low eigenvalues of the fermionic determinant. Fortunately the short distance behaviour of 
QCD is very well understood. This should allow the identification of the important degrees 
of freedom and lead to a method for including the effects of the higher eigenvalues into the 
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Figure 15: Quark-antiquark potential for quenched theory and dynamical theory 
with sea quark kappa values 0.1570, 0.1597. Error bars shown only for lightest mass 
(errors are largest for this case). 

Monte Carlo process. In particular, we know that for sufficiently high momentum scales 
this physics should be accurately described by an improved gauge action involving Wilson 
loops on short distance scales only. 

The fermion determinant can be separated into two pieces 

In detH = [Tr In H]io^ x + [Tr In H^.g^ a (20) 

where the lowest ricut eigenvalues are directly calculated and included in the Monte Carlo 
updating procedure. The contribution of the vast majority of the larger eigenvalues can be 
included by some approximation to the high end that (1) matches onto the low eigenvalue 
results without gaps or double counting, (2) is controlled and (3) becomes exact in the 
continuum limit. We can define the difference A between the approximate action, denoted 
Sa, and the exact contribution of the high eigenvalues, St = [Tr In H]high x as follows: 

A^St-Sa (21) 
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Any acceptable method must ensure that this difference is small (< 1) for each configuration. 
Therefore, we will demand that the variance of A for any set configurations is less than 
unity. Actually, we will find that a relatively simple effective loop action yields values of A 
considerably less than unity for interesting values of ricut- 

Two numerical methods suggest themselves for calculating the high eigenvalues of the 
fermion determinant: 

• The multiboson approach of Liischer|^. 

• Using a small number of gauge loops to model the determinant as proposed by Sexton 
and Weingarten jl^, and Irving and Sexton pO| . 

One method to compute the high eigenvalues which is guaranteed to succeed is the 
multiboson approach of LiischerQ. Define 

PeffiU) ^ [det(D + m)]"/ exp {-Sg{U)) (22) 

and 

= 75(D + m)/[c,n(8 + m)] (c„ > 1) (23) 

where Cm is chosen so that the eigenvalues of H are in the interval (—1,1). Consider two 
flavors of light Wilson-Dirac quarks {uf — 2). Liischer chooses a sequence of polynomials 
Pn{s) of even degree n such that 

lim P„(s) = 1/s for aU < s < 1 (24) 

n — ^oo 

then 

detff^ = lim [detP„(fl"2)]-i (25) 

Choose polynomials such that complex roots zi . . . Zn come in complex conjugate pairs (non 
real) so that ^/z ~ ^ ^ iu. Then 

n 

P(i/2) ^ const J] [{H ~fi, f + vl] (26) 
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and 

n 

detH'^= lim Y\ detUH - i^kf + ly?]-'^ (27) 

n— >oo 

fe=l 

Hence we can finally write 

Pe/ / (U) = lim ^ / D<PD<P^ exp -Sb{U, cj)) (28) 
where the bosonic action is given by 

n 

= Sg{U) + |(ir - i^k)4>k{x)\^ + A4>k{x)\\ (29) 

fe=l X 

To estimate how many boson fields are required to represent the original action to a fixed 

accuracy in the range (e < s < 1) Liischer considered polynomials of the Chcbyshev type 
(denoted T). Defining u = (s — e)/(l — e) and cos^ = 2u — 1, T*{u) = cos(to) and we can 
write 

P{8) = [l + pT*=i(«)]/s (30) 
R{s) = [P{s) - (l/s)]s(e < s < 1) 

where p is chosen so that the P(s) is finite as s — > 0. The error is given by: 

1^(^)1 ^2(i^)"+i (31) 

where n is the number of boson fields and the fit is cutoff for eigenvalues below e. Therefore, 
the convergence is exponential with rate as n ^ oo. 

The main practical problem with this multiboson method is that it requires an increas- 
ingly large number of boson fields as the quark mass becomes lighter. As iriq — > 0, we must 
take e — > 0, but to obtain a fixed level of accuracy we must hold 2,/en fixed and hence n 
increases without bound. 

However the multiboson method matches nicely onto the calculation of low eigenvalues 
discussed previously. In our truncated determinant method, the cutoff e for the multiboson 
method is set by the highest eigenvalue of s = (75 (-D + m))^ which is explicitly included in 
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the low end calculations. Hence it does not explode as the quark mass goes to zero. The 
combination of methods remain accurate for all quark masses. For example, for /3 = 5.9 on 
a 12^x24 lattice with direct inclusion of the lowest 100 eigenvalues, the associated cutoff for 
the multiboson simulation of the high eigenvalues is -v/e « 0.035 independent of the light 
quark mass. 

Furthermore, the error associated with the inaccurate behaviour of the polynomial fit in 
the range < s < e can be corrected as low eigenvalues are computed for every configuration 
update. We obtain a reweighting term. 



which can be included to eliminate errors in the region < s < e. 

Using the multiboson method for the high end of the determinant satisfies all our re- 
quirements and completes the algorithm. However it is interesting to study if we can reduce 
the total required computations even further using a more physical approach to the high 
eigenvalues. First, consider how many of the high eigenvalues we are computing actually 
have physical information and are not just lattice artifacts. For example, for a 12^x24 lat- 
tice with (3 = 5.9 and k = .1587 there are 497,664 total eigenvalues of the Wilson-Dirac 
operator. We can explicitly calculate the number of eigenvalues less than some high energy 
cutoff. Using 1 GeV we have approximately 1500 eigenvalues (0.3%). For a fixed volume V 
and quark mass niq only a decreasing fraction of the eigenvalues are below a fixed physical 
scale as /3 ^ 0. Therefore, most of the range of large s fit in Liischer's multiboson method 
is physically unimportant. 

This suggests a more physically motivated method for dealing with the high eigenvalue 
part of the fermion determinant, in which one approximates the ultraviolet contribution to 
the quark determinant with an effective gauge action: 




(32) 



i=l 




(33) 
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where each Li is a set of gauge Hnks which form a closed path. The natural expansion is 
in the number of links. For zero links we have Lq, which is just a constant, for four links 
we have a plaquette, and six links give the three terms found in considerations of improved 
gauge actions 

This idea was originally proposed by Sexton and Weingarten p9[ and studied in more 
detail by Irving and Sexton[Q. These studies were done on a 6'* lattice at /? = 5.7 with 
Hybrid Monte-Carlo full QCD simulations (with a heavy sea quark). Their results were 
rather discouraging. It was hard to get a good approximation to the determinant with a 
closed set of loops and they needed large loops to even approach a reasonable fit|pO|. 

There are however two important differences between their study and our situation. 

• They simulated the whole determinant, while here we only need to approximate the 
eigenvalues above some cutoff. Hence we would expect the small loops to dominate 
at least for sufficiently high cutoff. 

• They used an approximate procedure to estimate stochastically the logarithm of the 
determinant needed, while we are exactly computing all eigenvalues for this study. 

It turns out that these differences are critical, as using approximately the same lattices (and 
with even lighter quarks) we find an excellent approximation to the high end with only small 
loops. 

We generated a set of 75 configurations on a 6* lattice at (3 = 5.7 and k = .1685. We 
included the lowest 30 eigenvalues (which corresponds to a physical cutoff of approximately 
~350 Mev) in the Monte Carlo accept/reject step in the generation of these independent 
configurations. The spectrum of eigenvalues is shown in Fig. p^ . 

We can see the importance of the higher eigenvalues by separating the high and low part 
of the fermion determinant for each configuration. This is shown in Fig. ^ Unlike the case 
of QED2 (cf. Fig 1) it is apparent that the UV contribution to the determinant definitely 
involves large fluctuations. Of course, the issue here is just whether these fluctuations are 



29 




eigenvalue 

Figure 16: Spectrum of eigenvalues for 6^ lattice,/?=5.7,K=0.1685, Nx=30 

well described by a simple effective gauge action. 

Considering only the high eigenvalues, an excellent fit to the fluctuations is obtained 
including four and six link closed loops. The variance of the fit is 0.265. The comparison 
between the fluctuations in the exact and approximate actions for the high eigenvalue piece 
is shown in Fig 

As expected, if only the plaquette term is included the variance is larger (2.25) and we 
must move the low eigenvalue cutoff to = 50 (w 700 MeV) to reduce the variance below 
one. The results for various cutoffs and terms included are shown in Table 1. 

Although more study is required this second method looks very attractive for dealing 
with the high end of the fcrmion determinant in full QCD with light dynamical quarks. 
Simulations would be performed by including the predetermined effective gauge action Sa 
in the gauge updates and computing the infrared part of the determinant as in the truncated 
determinant simulations described in this paper. 

In summary, we have found that there are at least two viable methods to deal with the 
contribution of the fcrmion determinant not computed explicitly in a truncated determinant 
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Figure 17: Fluctuations in the low (n < 30) and high eigenvalues for the same 
configurations as previous figure. 
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Figure 18: Fluctuations in the high eigenvalue piece of Trlnff are indicated by a 
dashed line. Fluctuations in the fitted effective gauge action (6 links and less) are 
indicated with open circles 
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nx cut 


A (MeV) 


4 link fit 


6 link fit 


6 link fit 








wo WL 


with WL 








4.98 


1.074 


0.835 


±15 


340 


2.25 


0.2652 


0.233 


±50 


700 


0.940 


0.0564 


0.0491 


±250 


1,210 


0.0733 


0.0695 


0.0641 


±1250 


2,220 


0.138 


0.0198 


0.0180 



Table 1: Fits to the high eigenvalues of the quark determinant by various small 
gauge loops (WL denotes a Wilson line passing through the entire lattice) . 

approach. 

7 Conclusions 

We have proposed an algorithm for Monte Carlo simulation of full QCD with light dynamical 
quarks which has as its central feature a separation of the quark determinant into products 
over low and high eigenvalues. This separation is a direct reflection of the different physical 
roles played by these two sectors, with high eigenvalues (small quark loops) having the 
primary effect of modifying the gauge interaction strength, while the low eigenvalues (large 
quark loops) determine the long-range chiral structure of the theory. Our procedure for 
generating full QCD configurations then entails an exact calculation of some fixed number 
of the lowest lying eigenvalues oi H = 75(I/'(^)— m), using a Lanczos algorithm. This careful 
treatment of low Dirac eigenmodes is motivated by the conviction that the most essential 
differences between quenched and full QCD reside in their long-range chiral structure and 
associated topological properties. 

The complete formulation of our proposed algorithm also allows various methods of 
incorporating the higher eigenvalues which were omitted in the Lanczos calculation. For 
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example, the algorithm matches cleanly onto the multiboson approach. We have also ex- 
plored a particularly promising approach for incorporating the correct ultraviolet behavior 
by the use of a sum over relatively small Wilson loops to represent the high-eigenvalue con- 
tribution to InDet H. By doing a complete diagonalization on 6^ gauge lattices, we found 
that this Wilson loop representation works well with only a few small Wilson loops (4-link 
and 6-link) included. Here, the small loop approximation to the tnmcated InDet H only 
succeeds when the lowest eigenvalues (< 300-400 MeV) are excluded, so it forms an ideal 
complement to the Lanczos treatment of the low eigenvalues. In the Monte Carlo simula- 
tions discussed in this paper, we have used the pure Wilson plaquette gauge action for the 
heat-bath sweeps and carried out the accept-rcjcct test using the truncated low-eigenvalue 
contribution to the determinant. In the context of the more general Wilson-loop description 
of the high-eigenvalue segment, our present procedure is equivalent to approximating the 
high-eigenvalue contribution to InDet H by just the sum of a constant and a plaquette 
term. (This is in the same sense that the ordinary quenched approximation is equivalent to 
approximating the full determinant by a shift of /?.) 

Perhaps the best news to emerge from these numerical simulations is that the Metropolis 
test on the low-eigenvalue-truncated determinant yields a reasonably large acceptance ratio 
after one or more complete heat bath sweeps, even when the quarks are very light. This is 
in marked contrast to what would happen if one tried to include the full determinant via an 
accept-reject step between quenched Monte Carlo sweeps. Even for a single heat bath sweep, 
the fluctuations of the full determinant arc much too large to yield a reasonable acceptance 
rate. The Lanczos calculation of the lowest few hundred eigenvalues requires an amount 
of computing time of the same order as that of an ordinary conjugate-gradient inversion of 
the Dirac operator. Thus, even if the accept-reject step is performed after every sweep, the 
computing required is still comparable to that of other full QCD algorithms such as hybrid 
Monte Carlo. Moreover, the performance of the algorithm docs not seriously degrade in the 
light quark limit, which may provide a significant advantage over hybrid Monte Carlo for 
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the study of chiral behavior in hill QCD. Finally, for issues associated with chiral symmetry, 
the special handling of low eigenvalues is theoretically appropriate, and the eigenmodes 
extracted by the Lanczos analysis provide a detailed view of the connection between chiral 
symmetry breaking and the low-lying Dirac spectrum. 



34 



Acknowledgements 

A. Duncan is grateful for the hospitality of the Fermilab Theory Group, where this work 
was performed. The work of A. Duncan was supported in part by NSF grant PHY97-22097. 
The work of E. Eichten was performed at the Fermi National Accelerator Laboratory, which 
is operated by University Research Association, Inc., under contract DE-AC02-76CHO3000. 

The work of H. Thackcr was supported in part by the Department of Energy under grant DE- 
AS05-89ER40518. Much of the numerical work was performed on the Fermilab ACPMAPS 
system. 



35 



References 

[1] W. Bardeen, A. Duncan, E. Eichten, G. Hockney and H. Thacker Phys. Rev. (1998). 
[2] W. Bardeen, A. Duncan, E. Eichten, and H. Thacker, Phys. Rev. (1998). 
[3] W. Bardeen, A. Duncan, E. Eichten, and H. Thacker, |hep-lat 9806002| . 



[4] "The Uses of Instantons" , in The Aspects of Symmetry, Selected Erice Lectures of S. 
Coleman, (Cambridge 1985). 

[5] R. Setoodeh, C.T.H. Davies, and I.M. Barbour, Phys. Lett B213 195 (1988); K.M. 
Bitar, A.D. Kennedy, and P. Rossi, Phys. Lett. B234 333 (1990). 

[6] H. Leutwyler and A. Smilga, Phys. Rev. D46, 5607 (1992). 

[7] R. Narayanan and H. Neuberger, Nucl. Phys. B443 305 (1995). 

[8] A. Morel, J. Phys (paris) 48,(1987)1111; S. Sharpc, Phys. Rev. D41, (1990) 3233; C. 
Bernard and M. Golterman, Phys. Rev. D46 (1992)853. 

[9] M. Liischcr, Nucl. Phys. B418, 637 (1994). 

[10] J. Smit and J. Vink, Nucl. Phys. B286(1987)485; J. Vink, Nucl. Phys. B307(1988)549. 

[11] A. Duncan (speaker), W. Bardeen, E. Eichten, and H. Thacker, talk given at Lattice 
97 (Edinburgh), in Nucl. Phys. Proc. Suppl. 63, 811 (1998). 

[12] T. Kalkreuter and H. Simma, Comput.Phys.Commun. 93 33 (1996). 

[13] G.H. Golub and C.F. Loan, Matrix Computations, 2nd edition (Johns Hopkins, 1990). 

[14] J. CuUum and R.A. Willoughby, J. Comp. Phys. 44 329 (1981). 

[15] T. Kalkreuter, Comput.Phys.Commun. 95 1 (1996). 

[16] A. Duncan, E. Eichten, J. Flynn, B. Hill, G. Hockney, and H. Thacker, Phys. Rev. 
D51 5101 (1995). 



36 



[17] T. dcGrand and A. Hasenfratz, Phys. Rev. D49 466 (1994). 

[18] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in 
C, 2nd ed. (Cambridge University Press 1992). 

[19] J. C. Sexton and D. H. Weingarten, Phys. Rev. D55 4025 (1997). 

[20] A. C. Irving and J. C. Sexton, Phys. Rev. D55 5456 (1997); A.C. Irving, J.C. Sexton 



and E. CahiU, |hep-lat/9708004 , 
[21] M. Liischcr and P. Wcisz, Phys. Lett. 158B 250 (1985), and references therein. 



37 



